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SOLVING UPWIND-BIASED DISCRETIZATIONS II: MULTIGRID SOLVER USING 

SEMICOARSENING 


BORIS DISKIN' 


Abstract, This paper studies a novel multigrid approach to the solution for a second order upwind- 
biased discretization of the convection equation in two dimensions. This approach is based on semicoarsening 
and well balanced explicit correction terms added to coarse-grid operators to maintain on coarse grids the 
same cross-characteristic interaction as on the target (fine) grid. Colored relaxation schemes are used on 
all the levels allowing a very efficient parallel implementation. The results of the numerical tests can be 
summarized as follows: 

1. The residual asymptotic convergence rate of the proposed V(0, 2) multigrid cycle is about 3 per 
cycle. This convergence rate far surpasses the theoretical limit (4/3) predicted for standard multigrid 
algorithms using full coarsening. The reported efficiency does not deteriorate with increasing the 
cycle depth (number of levels) and/or refining the target-grid mesh spacing. 

2. The full multigrid algorithm (FMG) with two V(0, 2) cycles on the target grid and just one V 7 (0, 2) 
cycle on all the coarse grids always provides an approximate solution with the algebraic error less 
than the discretization error. Estimates of the total work in the FMG algorithm are ranged between 
18 and 30 minimal work units (depending on the target discretization). Thus, the overall efficiency of 
the FMG solver closely approaches (if does not achieve) the goal of the textbook multigrid efficiency. 

3. A novel approach to deriving a discrete solution approximating the true continuous solution with 
a relative accuracy given in advance is developed. An adaptive multigrid algorithm (AMA) using 
comparison of the solutions on two successive target grids to estimate the accuracy of the current 
target-grid solution is defined. A desired relative accuracy is accepted as an input parameter. The 
final target grid on which this accuracy can be achieved is chosen automatically in the solution 
process. The actual relative accuracy of the discrete solution approximation obtained by AMA is 
always better than the required accuracy; the computational complexity of the AMA algorithm 
is (nearly) optimal (comparable with the complexity of the FMG algorithm applied to solve the 
problem on the optimally spaced target grid). 

Key words, convection, upwind-biased discretization, multigrid solvers, textbook multigrid efficiency 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. Full multigrid (FMG) algorithms are known to be very efficient solvers for compli- 
cated systems of partial differential equations. It was rigorously proved (see [3] and [4]) that these algorithms 
can solve a general discretized elliptic problem to the discretization accuracy in a computational work which 
is only a small multiple of the operation count in the discrete problem itself. This efficiency is called the 
textbook multigrid efficiency (TME). The target-grid (grid h) solution obtained by a FMG solver is usually 
required to satisfy to the following condition: its algebraic error || u fl — u h |J must be less than the discretization 
error ||u /i -U h ||, where u h is the exact discrete solution, u h is the FMG solution, U h is a reasonable target- 
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grid representation of the true solution of the differentiable equation, and || • || is a given norm of interest. 
Then the total error \\u h -U h \\ is bounded above by sum of the algebraic and discretization errors. In solving 
nonelliptic problems, regular FMG algorithms sometimes fail to achieve such an accurate solution and then 
other very work-consuming methods are applied (either instead of or in addition to a FMG algorithm) to 
solve the problems. 

The goal of achieving TME in solving complicated computational fluid dynamics (CFD) problems at- 
tracts recently many researchers (see, e.g., the list of references in [5]). Many of textbook efficient multigrid 
CFD solvers developed in the last decade are based on the idea of separation of the elliptic and nonelliptic 
factors contributing to the corresponding system of partial differential equations. (See [9], [14], [15], [17], 
[18].) In all these solvers, the nonelliptic part was represented by the advection operator. 

The simplest way to solve the advection operator is to employ the downstream marching. If the cor- 
responding discretization is a stable upwind discretization and the field of velocities does not recirculate, 
then this marching proves to be a very efficient solver yielding an accurate solution to a nonlinear advection 
equation in just a few sweeps (a single downstream sweep provides the exact solution to a linearized prob- 
lem). However, if a discretization of the advection operator is not fully upwind (e.g., only upwind biased) 
the marching in its pure form is inapplicable. Other single-grid methods like defect-correction iterations or 
predictor-corrector technique can be used instead, in the recent paper [12], we showed that the number of 
defect- correct ion sweeps required to derive an accurate solution to the target discrete problem may be grid 
dependent. For example, in problems where the target operator is second order accurate while the correction 
is computed in solving a first order accurate discretization, the necessary number of iterations grows on fine 
grids approximately as ft -1 / 3 , where h is the grid meshsize. A detailed analysis of the predictor-corrector 
scheme is the subject of future studies. In many practical cases, these single-grid methods may result in 
very efficient solvers. However, each of them implements, one way or another, the idea of marching. In 
other words, their efficiency is significantly based on the correctness of the order in which values at grid 
nodes are updated. The consequent order of grid passages seems to be a serious obstacle for transferring 
these algorithms on parallel computers. In this paper, we suggest a novel multigrid solver to the advection 
equation employing colored relaxation schemes on all the levels. Coloring implies that discrete equations at 
grid nodes of the same color can be relaxed simultaneously (in parallel). Such schemes naturally possess a 
very high level of parallelism. 

It has long been known that standard multigrid solvers to the advection equation employing full coars- 
ening suffer an inherent slowdown in problems where the velocity direction does not coincide with the grid 
lines (see [1], [2]). This poor convergence is explained by an increased cross-characteristic interaction (e.g., 
dissipation) on coarse grids. The difficulty associated with the cross-characteristic interaction is very promi- 
nent in homogeneous problems with characteristics emanating from the boundary. In these problems, the 
quality of the coarse-grid approximation is determined by how well certain incoming oscillations are advected 
from the inflow boundary into the domain. The increased coarse-grid cross-characteristic interaction causes 
a decay and a phase shift of these incoming oscillations which are significantly different from their values 
on the fine grid. Following [6], the cure proposed in this paper is to use semicoarsening together with a 
well balanced correction of coarse-grid operators allowing the coarse grids to maintain essentially the same 
cross-characteristic interaction as on the target (fine) grid. The resulting multigrid cycles demonstrate good 
asymptotic convergence rates, far overcoming the theoretical limit for full-coarsening algorithms. 

The efficiency of the present multigrid algorithm does not deteriorate in multiple iterations, i.e. it 
demonstrates a good asymptotic convergence rate. Nevertheless, we believe that the role of a fast asymptotic 
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convergence is often overestimated. It is true that in elliptic problems the asymptotic convergence rate is 
defined by the slowest-to-converge error component and, therefore, a good asymptotic rate implies good 
convergence throughout the solution process. In solving nonelliptic problems, a fast asymptotic convergence 
is usually achieved only when all the troubling components have already been carried out of the domain of 
interest by many (skwer) iterations (see, e.g., [10] and [12]). Thus, fast converging iterations usually follow 
many those with much slower convergence rates. Moreover, in some problems, yielding a very accurate 
solution does not necessarily require a fast per-cycle convergence. These considerations have brought us 
to another important issue addressed in this paper which is criteria for stopping further computing when 
an accurate enough solution is achieved. A widely used criterion is sufficiently small residuals in the cor- 
responding discrete equations. To satisfy this criterion, the algorithm is generally expected to be iterated 
many times and a good asymptotic convergence is, thus, really important. However, it w T as realized in many 
experimental and theoretical w r orks that most of the computational time under this criterion is expended for 
approximating the discrete solution rather than the true solution of the differential problem. An approxi- 
mation to the differential solution within the discretization accuracy is reached long before the residuals in 
the discrete problem are reduced to a desired low level. Moreover, small residuals cannot, actually, ensure 
a good accuracy as w r ell. In nonelliptic problems, some slow'-to-converge error components are quite smooth 
along characteristics and possibly oscillating in the cross-characteristic direction, therefore, their residuals 
are very small comparing with other error components of the same amplitude. Thus, small residuals (say, in 
Lo norm) does not necessarily imply small amplitudes for these components. 

The well-knowm alternative to the small-residuals stopping criterion is the FMG approach, where the 
necessary number of iterations on each grid is defined in advance. This approach assumes an a priori choice 
of a target grid and a discretization on it, and, then, the FMG algorithm provides a target-grid solution 
approximation w r ith the algebraic error less than the discretization error. This is a good option, especially, 
when the accuracy of the target-grid discretization and convergence properties of the FMG algorithm for the 
problem of interest were previously established. In many practical cases, however, where the true solution 
is unknowm, either analytically or from an experience, engineers opt to drive residuals to the machine zero 
level in order to be sure that at, least the discrete problem is solved. Another difficulty associated with 
both the “regular FMG and ‘"small residuals 77 approaches is the need of an off-line analysis in order to 
establish the accuracy of the obtained solution. In many cases, this analysis indicates that either the desired 
accuracy has not been achieved because of a bad discretization error on the chosen grid and the problem 
should be resolved on a finer grid, or the desired accuracy could be reached on a coarser grid in much shorter 
time. In this paper, w r e propose another criterion indicating that a solution approximation possessing a 
desired relative accuracy has been obtained. This criterion is the comparison of solutions on different grids. 
Our experiments with an adaptive multigrid algorithm (AMA) using this criterion w r here the choice of an 
appropriate final discretization grid is an essential part of the solver are reported in Section 5 (see also [12]). 

In Section 2, we formulate the model problem and introduce the main ideas for the multigrid treatment 
of the cross-characteristic interaction. In the next Section 3, we define and test multigrid cycles which, then, 
are employed in framework of FMG solvers (Section 4) and AMA solvers (Section 5). 

2. Model Problem Description. 

2.1. Differential Equation. The model problem we study in this paper is the tw T o-dimensional (2D) 
constant-coefficient convection equation 
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(2.1) 


LU = (a-v)lJ = F(x,y), 

where a — (a\,a 2 ) is a given vector. The solution U(x<y) is a differentiable function defined on the unit 
square ( x,y ) 6 [0,1] x [0,1]. 

Let 0 be the nonalignment angle (another name which is common in CFD is the angle of attack), i.e., 
the angle between the vector a and the positive direction of the x axis; t = tan<^ = a 2 /a j is the nonalignment 
parameter. For simplicity, we assume a horizontal inclination a\ > a 2 > 0 and, therefore, 1 > t > 0. We call 
the x axis the reference axis. 

Equation (2.1) can be rewritten as 

(2.2) ^tf = /(*,y), 

where f(x,y) = F(x, y)/\a\, |a| = \/af + of, and £ = is a variable along the characteristic of (2.1). 

Equation (2.1) is subject to Dirichlet boundary conditions at the inflow boundary x = 0 and periodic 
conditions in the y direction. 


(2.3) U(0,y) = g(y), U(x,y) = U(x,y + 1), 

where g(y) is a given function. 

In the 2D constant-coefficient case which is studied in this paper, characteristics of (2.1) are straight 
lines defined by the velocity direction (characteristic lines). A function is called a characteristic component 
if it is much more smooth in the characteristic direction than in other directions. Possible extensions to the 
three dimensions (3D) and to variable velocity fields are briefly discussed at Section 6. 


2.2. Target-Grid Discretization. The approach we follow is to use a fixed Cartesian coordinate 
system independent of the characteristic direction. The problem (2.2)-(2.3) is discretized on the 2D Cartesian 
uniform grid with meshsize h in both the x and y directions. The target discretization grid is always assumed 
to be a uniform (square) grid. 

Let Ui y j. 2 be a discrete approximation to the solution U(x,y) at the point (x,y) = (i\h,i 2 h). To derive 
a proper discretization, we exploit the idea of a low- dimensional prototype introduced in [6]. Briefly, the 
low-dimensional prototype is a good discretization of the target (nonelliptic) differential problem on the 
grid induced by intersections of the Cartesian multidimensional discretization grid with the characteristic 
(low-dimensional) manifold (characteristic line in our case). (See Figure 2.1.) For our studies, we choose 
the low-dimensional prototype to be the (one-dimensional) second order accurate four-point discretization 
of the first derivative, corresponding to the Van Leer’s scheme with k = 0. 


(2*4) 4hy/l + & \ +3u*i,# 2 + Uij-2,i 2 -2f j — f i i ,*2 * 

The 2D discretization is obtained from the low-dimensional prototype by replacing function values at 
the ghost points (points with fractional vertical indexes) by weighted averages of the values at the vertically 
adjacent genuine grid points. The resulting narrow discretization is defined by 


L = ITT 


+f 2 /t 


(2.5) 


+ i .'*2 “b 3 Ui lf i 2 5u Z] _ i j 2 + tfq— 2,*2^ 

Ft(ui l+ i,j 2 +i 4- 3 Ui u i 2 - -f uq-2,,'2-2^ = 


u = l,2,...iV- 1, i 2 = 1,2, . . . N % N = l/h; 
wo,/ 2 = gib h), u- i, f2 = g'{i 2 h). 
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The outflow boundary conditions at i\ ~ N are discretized by the second order accurate narrow upwind 
scheme. 


(2.6) 


LhuN ** = 27T+^T 


^ _ (^ UN i 2 ~ 4wyv_i ,*2 + U N- 2,i a ) 

|^3ua^ 2 — 4ltyv — l,t' 2 — i -f Ujv- 2 ,i 2 -2 j ^ = fN,i 2 - 


In numerical experiments reported below, we used a corrected outflow discretization to provide the same 
cross-characteristic interaction as in the interior of the domain. See details in Sections 2.3 and 2.4. The 
discretization of the right-hand side function is / M ,* 2 = F(ijft,i 2 ft)/ \a\. Function g'(y) is an additional 
numerical boundary condition. In model problems, where the exact solution U(x,y) is known, one can 
define g'(y) = U{-h,y). 

The discrete scheme (2.5) is upwind biased, i.e., not a pure upstream scheme, in the interior since, for 
defining the operator value at the point (ii , * 2 ), the solution values at the downstream points (i\ + l,i 2 ) and 
(z'i + 1, i -2 + 1) are required. 


2.3. Cross-Characteristic Interaction. The cross-characteristic interaction introduced by a discrete 
operator can be quantitatively estimated as the coefficient of the lowest pure cross-characteristic derivative 
appearing in the first differential approximation (FDA) to the discrete operator (see [19]). The FDA to (2.5) 
taken for the characteristic component is given by 


(2.7) 


FD 


. 4 ( 1 ") = ^ 


- h 


. t(l-Q(l-2Q 

12VT+72 




.t(l - t)(l - 3t. + 3t 2 ) 

8Vl + t 2 


yyyy 


We included two (the second and the third order) terms into consideration in order to secure a more gener- 
ality. Usually, the second order term (related to d yyy ) determines the main part of the cross- characteristic 
interaction; at some slopes (t & 0.5), however, this term degenerates and the following third order term 
becomes important. Moreover, involving these two terms together makes possible to apply the same argu- 
ments to analyzing other (possibly higher order) advection schemes, e.g., Van Leer’s schemes corresponding 
to different k. 

The true measure of the cross-characteristic interaction should be calculated as a coefficient of a derivative 
with respect to the cross-characteristic variable 7 / = . In the constant-coefficient narrow-discretization 

case, however, this coefficient is just proportional to the factor of the corresponding y- derivative. The 
cross-characteristic interaction induced by the operator itself is referred to as inherent cross- characteristic 
interaction , to distinguish it from the explicit cross- characteristic interaction introduced below. 

Previous studies on different types of nonelliptic equations (see [1], [6], [9] and [11] ) have shown that 
the main difficulty in constructing an efficient multigrid solver is a poor coarse-grid approximation to the 
fine-grid characteristic error components. It was observed that a coarse-grid operator defined on a grid built 
by full coarsening (i.e., when all the coarse-grid meshsizes are twice as large as their counterparts of the 
fine grid) unavoidably introduces a too strong cross-characteristic interaction. On the other hand, a narrow 
discretization on a semicoarsened grid (only the reference axis meshsize is doubled) results in a coarse-grid 
cross-characteristic interaction which is weaker than required. However, we can supply this operator on 
the semicoarsened grid with additional (explicit) terms, so that the total coarse-grid cross- characteristic 
interaction would be exactly the same as on the fine grid. 



2.4. Coarse-Grid Discretization. The coarse grids used in this multigrid construction are rectangu- 
lar grids with fixed (integer) aspect ratios m = h x /h y , where h x and h y are the meshsizes in the x and y 
axes respectively. 

The preliminary (without explicit terms) narrow coarse-grid discretization on the grid with aspect ratio 
m is derived from the low-dimensional prototype (cf. (2.5) and (2.6)). 


(2.8) 


(2.9) 


& flx ' hy ^'U‘j u i 2 _ i'mhjy/l+t* ^ ,S ) ( w ii+l,* 2 +* + 3 w ij,* 2 ^ u i \ + ^0 -2,i 2 -2fc^ 

+ 6 '( w ii + l.* 2 +<*+n +3 u ii,h ~ 5' u U-l,h-(k+\) + — * 2 ,i 2 — * 2 ( Ar-hl ) ) ^ = fiuh'y 

L {h ' My) U\fj 2 = ^ (l " - s ) ^3wA/ ? f 2 — 4 ua/- 1,,‘ 2 _A- + UA/_2, /2 _2fc) 

+«^3wa/,<2 - 4«a/_ 1 ,*2-(A-+l) + w Af-2,i 2 -2(*+l))^ = /Af,i 2 ? 

/, = 1 . 2, . . . M — 1 , » 2 = 1,2,... A 7 , A/ = l/ft x = l/(m/i tf ), A r = l//i y ; 


where A: + s = mf, A- is integer, 0 < .s < 1. The first, differential approximation to the discretization (2.8) in 
the interior of the domain is 


(2.10) FDA(L {h * M ^ =dz- Hi 


a *(l-*)(l-2s ) 0 , j 3 S (1 — ' s )(l ” 3s + 3.s 2 ) 0 

Umy/lTt 1 yyy y SmVTTT* yyyy ' 


The first differential approximation to the outflow boundary condition discretization (2.9) is given by 

( 2 . 11 ) FDA(L {h " h » )y ) =df - Hi 


2 s{l-a)(l-2a) n , i3 .s(1-s)(1-3.9 + 3.s 2 ) 

* ^ SmvTTF + 4mVTT^ 


dy 


Comparing the coarse-grid FDAs (2.10) and (2.11) with the target grid FDA (2.7) one can derive the final 
form of the coarse-grid operators on a grid with an aspect ratio rn 


( 2 . 12 ) 


L^ lr ' h »)ui i j 2 — 4m A y ^/i+f2 ^ ( m *i + M2+* + 3 Uf lt i 2 bUii-i'iz-k + Wij-2,i2-2fc j 

+ *( u h + U2+(AH-l) +3 U iu i 2 ~ 5Wj 1 _ li .i a _( *+i) + W,: 1 _2,f 2 -2fA + l)^ 

^ 7T^ .*2+2 3Uij ,i 2 -f 1 + 3u j ll2 _ j — Uj li j 2 _2^ 

~^7rf .*2+2 — ,/ 2 -f 1 +6 — 4ui ly i 2 -i + Ujj ,z 2 — 2^ = fi \ ,i 2 7 

L (/|j "''“)«A/./ 2 - 2, n ft y b+/ 2 ^ (* “ S ) (^ u M,h _ 4um- 1,» 2 -* + f*Af-2,i 2 -2*) 

+ s(3»<M,t 2 — 4UA/_ lii2 _(A + i) + “M-2,ij-2(Hl))j 

d - f^ii ,i 2 +2 — 3'Ui lf i 2 + i + 3u, 1?2 _ 1 — Uj 1? j 2 _ 2 ^ 

T 7T^ a’ 2+2 ~ ^^u,*2 + l ~l"6u, l? j 2 — “b Uij ,i 2 _2^ = /A/,i 2 j 

ii = 1,2,... A/ - 1, 22 = 1,2,... A 7 , M = l/h x = 1 /{mhy), N = 1 /h y \ 
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O 2D discretization stencil 

• Low-dimensional prototype stencil 
Characteristic line 


Fig. 2.1. 2D coarse-grid discretizatton stencil 


where 

A:i — i 2v /r+^ ( — ^ 0(1 ~ 2f) + s(l — s)(l - 2 s)/rnj ; 

(2 

See Figure 2.1 for a pictorial representation of the discretization stencil in the interior of the domain. 

This choice of coarse-grid operator in combination with semicoarsening provides the same absolute pff 
and relative fiff approximation orders of to L h with respect to characteristic components, Ph = 

Ph = 2 (in terminology of [20]). This, in turn, ensures us a good coarse-grid correction to the characteristic 
error components. 

Discretization (2.12) implies a small correction to the target-grid discretization of outflow boundary 
conditions. Generally speaking, the cross-characteristic interaction in discrete operators at the target-grid 
outflow boundary does not play an essential role since its influence on the solution in the interior is minimal. 
On coarse grids, however, where outflow boundary operators are “responsible” for a larger part of the 
domain, it is important to adjust their total cross-characteristic interaction to that interaction in the target- 
grid interior operator. We decided to apply the outflow boundary discretization proposed in (2.12) already 
on the target uniform grid (rn — 1). It means that the target-grid outflow boundary discretizations include 
explicit terms with non-zero coefficients and A 4 . The values of and A 4 in the interior of the target 
grid are zeros by definition. 


A 4 = (t( 1 - t)(l - 3* + 3 t 2 ) - s(l - s){l - 3.s + 3 s 2 )/mj ; 

.43 = (~t( 1 - t)(l - 2*)/12 + s(l - s)(l - 2s)/(3m)); 

A 4 = -jj-f (t( 1 - f)(l - 3f + 3f 2 )/8 - s( 1 - «)(1 - 3s + 3 s 2 )/(4to)). 
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2.5. Strong Cross-Characteristic Coupling. Our multigrid construction employs semicoarsening 
and narrow coarse-grid discretization schemes supplied with explicit higher order terms (which are discrete 
approximations to the vertical derivatives of suitable orders) in order to maintain on the coarse grids the 
same cross-characteristic interaction as on the target (fine) grid. Then, the characteristic error components 
are well eliminated by the coarse-grid correction. The noncharacteristic error components must be reduced 
in relaxation. On several finest grids, where the direction of the strongest coupling approximately coincides 
with the characteristic direction, one can derive a pointwise relaxation scheme which reduces efficiently all 
the components oscillating in the characteristic direction. However, successive semicoarsening implies a 
fast decrease in the inherent cross- characteristic interaction on coarse grids and, hence, a fast increase in 
the weight of the explicit terms in the coarse-grid operators (since the total coarse-grid cross-characteristic 
interaction remains fixed). Thus, the direction of strongest coupling tends to be vertical. After several 
semicoarsening steps, hence, any pointwise relaxation scheme fails to reduce the noncharacteristic error 
components. 

Roughly speaking, pointwise relaxation schemes loose their smoothing properties on grids where the 
“cross-character istic 71 coupling becomes stronger than the “characteristic 17 one. The qualitative description 
of relations between the couplings can be derived from the analysis of the symbol of the corresponding 
coarse-grid operator (see Section 3.3.1). This analysis tells us that the efficiency of pointwise relaxations 
may degrade on grids with aspect ratios m > 54. 

There are several ways to prevent this degradation: 

1. The first is to use vertical line relaxations in which points located on the same vertical grid line are 
relaxed simultaneously. Such a relaxation is required only on coarse grids with aspect ratios rn > 54. 
Practically it means that on fine grids (m < 32) a pointwise relaxation scheme can be efficiently 
used, while on coarser grids (m > 64) a line relaxation is employed. This method can be efficiently 
extended to the three dimensions (3D) (with a corresponding replacement of line relaxations with 
plane relaxations) and variable coefficient problems. This first method is extensively studied in this 
paper. 

2. The second method is to widen gradually the basic narrow coarse-grid discretizations in order to 
increase the inherent cross-characteristic interaction and to reduce, in this way, the weight of the 
explicit terms in the coarse-grid discretization. This method is very flexible. It eliminates the need 
of using line relaxations on coarse grids and, therefore, it is well suited for multiblock grids required 
for calculations in complex geometries. This approach is briefly discussed in Section 6. 

3. The third way is to use conditional coarsening technique, where a strong cross-characteristic coupling 
can be avoided by replacing part of the semicoarsening steps bv full coarsening steps (see [11]). This 
conditional coarsening approach seems to be slightly cheaper than other in computing time but 
considerably more complicated to program, especially in extensions to variable coefficients. 

3. Multigrid Cycles. 

3.1. Multigrid Cycle for Low-Dimensional Prototype. The first necessary step is to derive an 
efficient solver for the prototype problem. This solver, then, wall serve us as a model for a multigrid solver in 
the full dimension. The one-dimensional multilevel V(y\ , ia>) cycle we have studied for the prototype problem 
consists of a colored relaxation scheme, an upwind-biased residual transfer and a linear interpolation of the 
coarse-grid correction. On each level, except the coarsest one where the problem is directly solved, v\ 
relaxation sweeps are performed before transferring residuals to the coarse grid and i / 2 sweeps are performed 
after receiving coarse-grid corrections. Below , we present a detailed description of all the components of this 


8 



cycle. 

On the grid induced on the characteristic line, the one-dimensional prototype problem can be rewritten 
in new variables (cf. (2.4) and (2.5)) as 


(3.1) 


L {1D ^iiii = + 3 Ui x — 5n. ?1 _i — /m , 

*i = 1,2,. ..TV- 1; 

L (ID) uk = ^ \3u N - 4?t lV _] + u.v-i'j = In; 


u 0 = go, u-i = g i; 


where he is the meshsize between nodes of the low-dimensional prototype discretization. 

3.1.1. Relaxation Scheme. The downstream pointwiso Gauss-Seidel relaxation scheme for the one- 
dimensional prototype interior equation (3.1) is unstable. In spite of the fact that an immediate error 
explosion can be prevented by using a color order of the relaxation passage, we decided to pick a stable 
relaxation scheme. This is done mainly because of the observation that (high-frequency) instabilities in 
schemes are usually accompanied with bad smoothing properties (especially, in further multidimensional 
extension). The relaxation scheme chosen for the problem (3. 1) is a three-stage defect-correction-type scheme. 
Stage 1: The residual function for the target problem (3.1) is computed. 


(3.2) 


(i = fh ~L (1D) u h , 


ii = 1,2,... A' - 1; 


where is the current solution approximation. 

Stage 2: The correction i’,, is calculated by relaxing the system 


LdVi i =r u , ij = 1,2, ... TV - 1; 


where vj t (i\ = - 1 , 0, . . . , TV, N + 1 ) are initialized by zeros and operator L d is the target, operator 
C 1 r>> supplied with an additional (third order) dissipation term: 


LdVh 



(v il+ i + Zv n - 5t’ti — 
Uii+2 — 4w; l+ i +6u„ 


i + 


v u 


4i'j, -i + . 


The parameter A is chosen to maintain stability in the downstream marching. In our tests, we have 
picked A (approximately) equaled to the minimum positive value ensuring that all the (possibly 
complex) roots of the quadratic equation 


(~ + 6A)x 2 — ( — + 4A):r + (- + A) — 0 
4 4 4 

are inside the unit circle. It was found numerically, that A « 0.084. 

Stage 3: The current approximation Ui 1 is corrected to the improved approximation u u by 


Uii + v tl , i i = 1,2,. 

UN = ^ ^2/jv/^ - (-4«Af-l + 


..N- 1; 


We still have some freedom in choosing the order of the relaxation passage on Stage 2. It is proved 
that the best smoothing is observed in the downstream marching. However, smoothing rates of colored 
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relaxation schemes are also quite good. A colored (with p colors) relaxation order is defined as following. 
One colored relaxation sweep consists of p passages, each passes through (approximately) N/p points. In 
the first passage, all the points with coordinates v 1 = 1 + jp (j is a non-negative integer) are relaxed; in the 
second passage, all the points with coordinates t'i = 2 + jp are relaxed (in this passage the new values at 
previously relaxed points are used); and so on till all the points are updated. The minimum p required to 
enable a full parallelization and preclude the appearance of relaxation boundary layers is p = 4. This choice 
of p already provides a very good smoothing rate (see the 2D mode analysis in Section 3.3.2 where this ID 
prototype appears as a particular case of alignment). Moreover, the soothing rates is further improved for 
P > 4. Thus, in practical problems the maximal possible p should be picked on. In our tests, we have mostly 
experimented with p = 4. 

3.1.2. Residual Transfer. At this stage, we compute a coarse-grid approximation to the current 
fine-grid residual function (3.2). This fine-to-coarse transfer (restriction) is defined by 

(3-3) R it = i(r 2l , +r 2l| _,), 

where R and r denote the coarse- and fine-grid residual functions respectively. 

.. The coarse-grid correction V is interpolated (pro- 


Vi„ 

^(lq-1 -t r„). 

where v is the correction to the tine-grid solution approximation. 


3.1.3. Coarse-Grid Correction Interpolate 

longatcd) to the fine grid by 


(3-4) 







3.1.4. Efficiency. We have experimented with a two-level V(0,3) algorithm on different grids with 
different right-hand side functions /. The tests have demonstrated good convergence rates. The worst rate 
observed in a cycle during the solution process is about 3, while the asymptotic convergence rate is better 
than 6 per cycle. For p = 8 and the corresponding values are 3.7 (the worst per-cycle convergence rate) 
and 10 (the asymptotic convergence rate). The two-level discrete half-space and matrix analyses which are 
too cumbersome to be presented here (see some examples of the analyses in [8], [10], and [12]) predict the 
low bound for the convergence rates to be 2.3 per cycle. This bound does not depend on h and p providing 
p < h ~ 1 . 

3.2. Two-Level Cycles. In this section, we discuss the basic parts of the full-dimensional (2D) multi- 
grid cycle such as relaxation schemes, residual transfer and coarse-grid correction interpolation. 

3.2.1. Pointwise and Line Relaxation Schemes. The relaxation schemes defined in this section are 
derived from the colored one-dimensional schemes described in Section 3.1.1. The operators involved are, of 
course, extended to the two dimensions. The target operator L {h ^ h ^ on an anisotropic grid with an aspect 
ratio in is given in (2.12). The “driver” operator L d which is relaxed for the error equation on Stage 2 (see 
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Section 3.1.1), is defined as 


j ( h X Jly ) 

u d 




+ 



L( h *' h y)v iil2 

^M+2,i 2+2A* 4t^] + 1 ,<2+^’ d" _ i f t2 — k "h —2, #2 — '2k ^ 

v *i+2,* 2 +2(Jt+l) - ^ u i 1 + 1 ,i 2 + (*+ 1 ) 

+6t’t , ,i 2 — 4?;^ - 1 ,i 2 - ( k+ 1 ) + t'i , -2,12 -2( fc+1 ) ) 

(l - 3.s + 3s 2 ^ ^V / 1( * 2 +2 - 4n lt? ,; 2+ i + 6<;* 1t?2 - 4^.*, ,* 2 — 1 + t^ 1 ,i 2 


» 


where the nonalignment parameters k and s are defined in Section 2.4, A = 0.084, h $ ~ mh y y/ 1 + 1 2 , and 
the added artificial dissipation term approximates the differential operator A/f|c^^. 

The number of colors p determines the order of relaxing the vertical grid lines. The lines with horizontal 
coordinates i\ = 1 + jp, ( j € Z) are relaxed first, then the lines with i \ = 2 -f jp< and so on; the lines to be 
relaxed last an 1 those with ?‘i = (p — 1) -f jp. In all the numerical tests below, we used p — 4. The difference 1 
between line and pointwise relaxation schemes to be used is how the solution values at the same* vertical 
grid line are updated. In the line relaxation scheme, all the equations centered at the same vertical grid line 
are solved all together. Simultaneous replacement of solution values at all the grid nodes belonging to the 
line reduces residuals on this line to (nearly) zero. In a pointwise relaxation, the solution approximation is 
changed in a point to satisfy the only discrete equation defined at this point. The orders of relaxing the 
grid points on a line can be different. We use the four-color order in which relaxation on a line is performed 
in four passages. Each passage updates every fourth points. The first passage starts from the point with 
vertical coordinate i 2 = 1; the second, from the point with / 2 = 3; the third, from the point with t 2 = 2; 
and the last, from the point with ?* 2 ~ 4. This pointwise (sixteen-color) relaxation scheme is efficient and 
especially attractive for further implementations on parallel computers. 


3.2.2. Intergrid Transfers. The residual transfer to the semicoarsened grid is given by 
(3-5) Ri j „ j 2 = -5 ^2*, ,* 2 T (1 — s)r 2 *,_ ij. 2 ~h + sr 2l ; 1 -i ? i 2 _(^ + i)^ * 

where r, 1</2 = — V hx ' hy) u lu i 2 is the fine-grid residual function, and R-i lf i 2 is the coarse-grid residual 

function. This restriction operator possesses the low-frequency order fhp — 1 and the high-frequency order 
nift = 1 (see definitions of intergrid transfer orders in [20] and nole that in semicoarsening algorithm only 
two fine-grid component are coupled on the coarse grid). 


The coarse-grid correction operator is a linear interpolation defined by 

(3-6) ( ^‘•' 2 = *V 2 ’ . , 

[ V2ii-Ma — 2 + *1 ?1 - 1 ,i 2 -(*+! ) + U — 5 )^ti,i 2 +/k + sV/, ,,' 2 + (Ar+l ) J ? 

where V is the coarse-grid solution and v is the correction to the fine-grid solution approximation. The 
low-frequency and the high-frequency orders of this prolongation operator are rhp = rn p = 2. 

Thus, the intergrid transfers satisfy the necessary conditions (derived in [20] and earlier in [4] and [13]) to 
provide a grid-independent convergence: nip + nip > 1 and nmx(inp, nip) > 0. 

3.2.3. Numerical Tests. In the full dimension, a two-level cycle V 2 (iq,i/ 2 ) employing semicoarsening 
can be defined as the following six steps 
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Step 1 Prerelaxation sweeps. The current approximation is improved by v\ relaxation sweeps. 

Step 2 Residual transfer . The coarse-grid approximation to the fine-grid residual function is calculated by 
means of (3.5). 

Step 3 Coarse- grid operator. The new values of the coarse-grid parameters rn, k . «s, A :1 , A 4 , A3, and A 4 are 
calculated. The coarse-grid equations (2.12) are formed. 

Step 4 Coarse-grid solution . The coarse-grid problem is solved. Its solution approximates the fine-grid 
error function. On this stage, we do not specify the solution method. It can be any method (direct 
or iterative) allowing to obtain an accurate solution to the coarse-grid problem. 

Step 5 Coarse-grid correction. The coarse-grid solution is interpolated by (3.6) to the fine grid. The current 
fine-grid approximation is corrected. 

Step 6 Postrelaxation sweeps. The current fine-grid approximation is improved by v 2 relaxation sweeps. 

We restricted ourselves to considering l 2 (0,2) cycles. Such cycles (with - 0) are especially attractive 
for parallel computing (see [7]). In numerous computational tests, we compared the performance of l 2 (0, 2) 
cycles with either full or semicoarsening on grids with different aspect ratios m = 1, 2, 4 , . . . , 512, 1024. 
Within the cycles employing semicoarsening, two types of relaxations (pointwise and line relaxations) were 
tested. 

In all the two-level tests, we used the zero right-hand side function f u j 2 . The inflow boundary conditions 
were chosen so that the function U (:r, y) = sin(u ;{y - tx)) is the exact continuous solution of the homogeneous 
problem (2.2), (2.3). The coefficients A3, A 4 , A 3 , and A 4 of the explicit terms in the fine-grid discretization 
were derived from the assumption that this fine grid itself was obtained from a uniform grid by (log 2 m steps 
of) semicoarsening. In other words, the total cross-characteristic interaction in the fine-grid discretization 
was the same as the inherent cross-characteristic interaction in the interior of a uniform grid with meshsize 
h y . 

A representative sample of the experimental results is shown on Figure 3.1. Each experiment included 
three different runs, each starting from the same initial approximation obtained by interpolation from the 
solution on semicoarsened grid. Both the pointwise and the line relaxation schemes used in the runs are 
described in Section 3.2.1 above. Run 1 (marked by pluses) employed the pointwise relaxation scheme and 
full coarsening. For small aspect ratios, where the inherent coarse-grid cross-characteristic interaction was 
stronger than the desired total cross-characteristic interaction, the coarse-grid values of A3, A 4 , A3, and A 4 
were set to zero. Run 2 tested the pointwise relaxation scheme together with semicoarsening. Run 3 used 
the line relaxation scheme and semicoarsened coarse grid. Each run consisted of 20 cycles; the ratio of 
norms of the residual before and after each cycle is calculated and pictured in the diagrams on Figure 3.1. In 
all the graphs on Figure 3.1, the horizontal coordinates serve to mark the cycle numbers. On vertical axes, 
the per-cycle convergence rates are displayed. Other notations used in titles on Figure 3.1 are following: 
rn — h T /hy is the fine-grid aspect ratio, where h x and h y are the corresponding x- and ^/-directional fine-grid 
meshsizes; t is the nonalignment parameter; uj is the frequency of the incoming component; RC is a relative 
coupling parameter. RC > 1 corresponds to discretizations for which smoothing properties of pointwise 
relaxations deteriorate. The methodology of calculation this RC parameter is described in Section 3.3.1. In 
all the tests pictured in Figure 3.1, we picked the nonalignment parameter t ~ 0.2, except the last test where 
t = 0.98. The value t = 0.2 roughly corresponds to the maximum inherent cross-characteristic interaction. 
The frequencies u) were chosen to satisfy to the two conditions: first, u> ~ 2 n 7r, (n is integer) to reduce the 
total computational time exploiting the vertical periodicity; second, the fine-grid discrete solution should 
provide a reasonable accuracy in approximating the true solution of the differential equation. The last 
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= 2 “ , 


32 




h y = 2 15 , m = h x /h y = 512, t = 0.2, RC= 15.925, co= 1024ji. 



Fig. 3.1. Residual convergence history m two-level experiments: Run i (pluses) is Vb(0, 2) cycle with full coarsening; 
Run 2 (triangles) is 12(0,2) cycle with semicoarsening and pointwise relaxation; Run 3 (pentagrams) is \'2(0,2) cycle with 
semicoarsening and line relaxation. 


experiment was added to confirm the claim that the algorithm efficiency, actually, depends on the relative 
coupling RC rather than on aspect ratios. (Compare the last experiment and the test on grid with m = 32). 
In particular, in cases of alignment (£ = 0or£ = l), the pointwise relaxation scheme can he efficiently applied 
on any grid with any aspect ratio. 

The first very prominent observation from analyzing the results of numerical tests is the superiority of 
the semicoarsening algorithms over the algorithm with full coarsening. The two semicoarsening algorithms 
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show a similar behavior for discretizations with small values of the relative coupling (RC < 1, m < 32). 
When RC > 1, only the algorithm employing the line relaxation scheme demonstrates a fast convergence. 
Note 1 Attempts to use higher order intergrid transfer operators result in some improvements in the full- 
coarsening algorithm, especially, on grids with RC « 1. However, this improved convergence still 
cannot compete with convergence rates demonstrated in the semicoarsening algorithm employing 
the line relaxation scheme. 

Note 2 Inclusion of these higher order operators into semicoarsening algorithms deteriorates the convergence. 

3.3. Fourier Mode Analysis. 

3.3.1. Coupling Analysis. The coupling analysis presented in this section studies the discrete opera- 
tor symbol L(0), 0 = (9 x ,9 y ). By definition, the symbol of a discrete operator is the response of this operator 
on the discrete Fourier component 


J^{h x J\ y) e i\9 x i\+9 y i- 2 ) _ £^ e i(0xii+0 y i2) 

The smoothing properties of pointwise relaxation schemes for discretized partial differential equations are 
essentially determined by the measure of /i-ellipticity in the target discretization. (See [1] and [2].) Briefly, 
the measure p of /i-ellipticity is calculated as i = min \L(9)\, where the minimum is taken over all the high- 
frequency Fourier components. The Fourier mode with normalized frequencies (|tf x | < 7 T, \0 x \ < i r) 

is called a high-frequency mode if max(|® x |, |0*|) > tc/2. A discrete operator is called /i-elliptic, if // is 
separated from zero (// > constant > 0). Roughly speaking, the larger absolute values of the discrete 
operator symbol for given high-frequency Fourier components the better these components are eliminated 
by relaxation from the error function. 

The multigrid construction proposed in Section 3.2 ensures a good approximation to the characteristic 
components of the solution. The noncharacteristic error components must be removed in relaxation. Let 9^ 
be a normalized characteristic frequency defined by 0$ + x = (9 X + mtO y ) mod 27r. The noncharacteristic 
components correspond to |0$| > 7r/2. Thus, while the discretization is semi h-elliptic in the characteristic 
direction, i.e., \L(9)\ is large enough for all 0 satisfying \9^\ > 7r/2, one can derive a pointwise relaxation 
scheme which efficiently reduces the noncharacteristic error components. The measure of the characteristic 
coupling (see Section 2.5) can be computed as 

fio = min \L(9)\. 
l * d >*/2 

In our constant-coefficient model problem, the cross- characteristic coupling is mainly determined by the 
vertical interactions and, therefore, its measure can be defined as 

u\ — min \L(9)\. 

\ 0 y \>*/'2 

On grids with large aspect ratios m > 64 where p 0 , the discretization becomes semi /^-elliptic in the 

vertical rather than in the characteristic direction and pointwise relaxation smoothing factors deteriorate 
for the noncharacteristic error components. Thus, the range of applicability of pointwise relaxations can be 
described by the value of the relative coupling parameter RC = pi/po: if RC < 1, then one can derive an 
efficient pointwise smoother, otherwise, a line relaxation should be used. The RC values in different two-level 
tests are displayed on Figure 3.1. 
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3.3.2. Smoothing Factor of Four-Color Line Relaxation. The symbol of the four-color line relax- 
ation Z(0) defined in Section 3.2.1 is a 4-by-4 error amplification matrix acting on a four of Fourier inodes 
e r(0iri+0yi 2 ) , j — o, 1,2,3; where 0{ are normalized frequencies satisfying 6{ + n = (0 X + jV/4) mod 2n 
(l&rl < 7T)* 


Co(8 w ) Ci(8 w ) 
C 3 (0 (I) ) Ca(B {l) ) 
C 2 {0 < 2 >) C 3 (0^) 
\ Ci(8 (3) ) C 2 (0< 3 <) 


C 2 (0 m ) C 3 (0 io) ) \ 
Ci(tf (,) ) C 2 (0 (1) ) 
Co(8 V2) ) Ci(0 (2) ) 
C :i (0< 3 >) C o (0 {3) ) / 


where — (0 J x ,O y ). Parameters Cq(0), C\ (0), C’2(0), and Cs(0) are defined as 


cm \ 



( 1 ) 


f V o (0) ) 


cm 

= M 


i 


Vi(9) 


cm 



i 


vm 


cm J 



1 J 


\ vm ) 



M is a constant, matrix 


i i i i \ 

-i —1 i 1 

-11-11 

i —1 —i 1 j 


Vo(9) = L(0)/L 3 (0 y ), 

v, (0) = ( L(0) - e- ie *L 2 (6 y )W0))/L 3 (e y ), 

V 2 (6) = ( L(9) - (0 y )V o (0) - e~ ie * L- 2 (9 y )\\ (S) 

-e™*L 3 (6 y )\o(8))/L 3 (0 y ), 

V 3 (8) = ( L(B) - e~ 2Wx Li(0 y )l\(9) - er' 9 * L 2 (0 y )V 2 {6) 
-e ie *L A (9 y )V o (0) - e 2ie * L 5 (6y)\'i (0)) /L 3 (0 y ); 
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L\ (0 y ) = I 1/4((1 - s)e~ i2ke « +ae- n (* + i) ff ») 


\ 

+a((1 — s)e~ l ' 2ke « + .se _t2 (* +1)e ») J /(my/T+l 2 ), 

L->(8 y ) = ^ — 5/4((l — s)e~‘ key + 

-4A((1 -8)e~ iM ' +se- i < H1 )"»)j/(mv/r+?), 

L:i(8 y ) = ^ 3/4 + A3i(sin(29 y ) — 2 sin(6^ y + .1( (2 cos(2^ v ) — 8cos(6 y ) + 6 ) 

+a(g - s(l - .s)(l - 3.s + 3s 2 )(2cos(20 ! ,) - 8cos(0„) + G)) j /(m\/TT7 2 ), 
L,,(0y) = ^ 1/4 ((1 - s)e ike » + se ,{k+}>6 ^ 

-4A((1 - s)e ik9 » + ««.*(*+> )«») J /(mVl + t 2 ), 

L T ,{8 y ) = ( a((1 — s)e r2kt>y + se ,2(fc+1|6, !'^ N j /{m\J 1 + t‘ l ). 


t is the nonaligiiment parameter. .4 3 , .4 4 , A, A\ and s are the parameters of discretization (2.12) on the grid 
with aspect ratio m. L(0) is the discretization symbol: 

L(0) = e-W’L^+e-iO'LMy) + L 3 (0 y ) + e' e *L A (6 y ) + e™* L 5 (9 y ). 

Following [16] we define the smoothing factor of the four-color line relaxation as the spectral radius of 
the matrix product Q(0)Z{9), where 


Q(9) = 


q 0 0 0 0 

0 q x 0 0 

0 0 q 2 0 

0 0 0 q 3 


qj = 1 U = !>• C 2, 3), if tt/2 < |0< j) | < t r, where (9‘ 2) + 7 r = ( 0 ( x j) + mt.6 „) mod 27r; 
and qj = 0 if |^ 2l | < 7r/2. 

We picked A = 0.084 and calculated smoothing factor Srn \ for a variety of different slopes t and aspect 
ratios m (all other parameters are derived from these). In all cases, Srrii < 0.54. The smoothing rate of 
two successive four-color line relaxation sweeps Sm *> which is defined as the spectral radius of Q(9)Z(0 ) 2 is 
Sm -2 < 0.45. 

3.4. Multilevel Cycles. We performed many experiments with a multilevel V r (0,2) cycle on uniform 
grids varying the nonalignment parameter t, the frequency uj of the incoming Fourier component, and the 
right-hand side function /. The multilevel cycle Vrf(i/i,i/ 2 ), where d is the cycle depth is defined similar to 
the two-level cycle (see Section 3.2.3) but Step 4 is replaced with recursive call of the same cycle applied to 
the coarse-grid problem. The coarsest grid where the problem is solved precisely is always the grid having 
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Fic;. 3.2. Residual convergence history in multilevel experiments on target uniform (square) grids. 


eight mesh spaces in the x direction in the interior part of the domain. We simplified the criterion of 
switching from pointwise to line relaxation to the following rule: pointwise relaxation sweeps are employed 
on grids with aspect ratios rn < 32; on grids with higher aspect ratios the line relaxation is used. Figure 3.2 
demonstrates results of numerical tests performed on different uniform target grids. Similar to the two-level 
tests we considered the homogeneous problem with boundary conditions derived from the assumption that 
function b (x, y) = sin (uj(y — tx)) is the exact solution of the differential problem. Initial approximations were 
obtained by interpolation of the exact discrete solution from the next immediate coarse grid with aspect ratio 
m = 2. In all this experiments, we picked t = 0.2 and solved the problems for incoming components with the 
same normalized frequency uh = tt/ 8. The graphical results confirm that the multigrid cycle convergence 
rates are grid independent. The residual convergence rate histories on the three finest of the tested grids are 
practically undistinguished. 

We also performed many numerical experiments testing the dependencies on the nonalignment parameter 
t and incoming frequencies u). The results can be summarized as following: 

1. In the first cycle, the L 0 0 residual norm of the error is reduced by almost an order of magnitude. 

2. In the first three cycles, the overall error reduction (also in the L 0 0 residual norm) is more than two 
orders of magnitude. 

3. The convergence rates in further cycles are ranged between 2 and 3.6 per cycle. 

4. Full Multigrid Algorithm (FMG). The goal of an FMG algorithm can be formulated as fast 

obtaining an accurate solution for a given discretization on a given target grid. A solution approximation is 
considered to be accurate if its algebraic error is less than the discretization error. In this section, we present 
an FMG algorithm based on the multilevel l 4 * * 7 (0,2) cycle. The setup work required for this algorithm can be 
described as the following four steps: 

Step 1. Target-grid problem. The discrete problem (2.12) is formulated on a uniform (m = 1) target grid. 

The total cross-characteristic interaction for the entire algorithm is defined as the inherent cross- 
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characteristic interaction in the discretization in the interior of the target grid. Proper discretizations 
of the right-hand side function and the boundary condition functions are also performed. In our 
implementation, the corresponding discrete functions /,, . g l2 and g[ 2 are formed by injection from 

their continuous counterparts. 

2. Next coarse grid. The next coarse grid is constructed by semicoarsening, where only meshsize in the 
reference x direction is doubled. 

3. Coarse- grid problem. The coarse-grid right-hand side function f c is formed bv the following aver- 
aging operator: 

■fU 

f M C ,*2 
*1 = 1 , 2 , 

where M ( and AI (M = 2M C ) are parameters defining the number of mesh spaces in the x direction 
in the interior of the domain on the coarse and fine grids respectively. N is the number of mesh 
spaces in the y direction (the periodicity direction). In semicoarsening, N is the same on all the 
grids involved in calculations. The parameters k and s are the fine-grid nonalignment parameters 
(k + s = mt, k is integer, 0 < s < 1, m is the fine-grid aspect ratio, t = tan <j) is the tangent of the 
nonalignment angle). The coarse-grid inflow boundary conditions are injected from the known true 
solution to the continuous problem. (See Note 1 below.) New values of coarse-grid discretization 
parameters such as the aspect ratio rn c , the nonalignment parameters k c and s c ' , the coefficients 
of the explicit cross-characteristic interaction terms in discretizations in the interior (A% and .4^ ) 
and at the outflow boundary (A% and ,4 £’) are calculated. 

Step 4 - Recursion. The Steps 2 and 3 are repeated until the coarsest possible grid is reached and its problem 
is defined. 

Note 1: Injecting true continuous solution values into the coarse-grid inflow boundary discretization is 
the easiest way to separate the issue of developing an efficient multigrid solver from the issue of deriving a 
proper high-order discretization to inflow boundary conditions. One possible solution for the latter is to use 
a central (or a downwind) discretization at the grid nodes adjacent to the inflow boundary. In this case, the 
coarse-grid inflow boundary conditions can be derived from those on the fine grid either by injection or by 
averaging. 

The execution of the FMG algorithm involves the following four steps: 

Step 1. Coarsest-grid solution. The problem on the coarsest grid is solved by some (direct or iterative) 
method. 

Step 2. Initial fine-grid solution approximation. The initial approximation u on the currently fine grid is 
derived from the coarse-grid solution u c by an interpolation which is the fourth order in the interior 
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of the domain (and the second order near the outflow boundary) in the characteristic direction. 
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Step 3. V Cycle. The obtained initial approximation is improved by one l r (0, 2) cycle. 

Step 4 . Recursion. Steps 2 and 3 are repeated until the target (finest) grid is reached. There one additional 
V(0. 2) cycle is performed. 

The total cost of this algorithm is about 30 minimal work units, where minimal work unit is defined 
as the number of computer operations required to evaluate residuals on the target grid. This work-unit 
count is about five times larger than usual in uniformly elliptic problems. It is contributed bv (1) somewhat 
expensive relaxation schemes using defect-correction type iterations, (2) semicoarsening increasing the cost 
of coarse-grid calculations, and (3) the need to perform the second target-grid cycle. In fact, if the target 
discretization was changed to the discretization used in the correction step within relaxation (which has the 
same approximation order as our target discretization), then the cost of relaxation sweeps would be twice 
as cheap and, therefore, the total of the whole FMG algorithm would be reduced to 18 minimal work units. 

In our numerical tests, we chose the right-hand side function / and and the inflow boundary condition 
function g so that the function U(x,y) — sin(^ J .x + 0 v y) was the exact solution of the continuous problem 
(2.2), (2.3). The tests were performed with a six-level FMG algorithm solving the problem on a uniform 
target grid with h x = h y = h = 2“ 8 . We experimented with different values of parameters 0 X and 0 y . For 
each component, we checked five different characteristic inclinations t = tan cj). Some representative results* 
are collected in Tables 4.1 and 4.2 where? the target-grid discretization error is compared with the algebraic 
errors at three stages: immediately after obtaining the initial target-grid solution approximation and at the 
end of the first and the second improving target-grid cycles. In the tables, i\ is the characteristic frequency 
0i t9y , and is the characteristic meshsize h$ = vTT t 2 h. A small absolute value of the normalized 

characteristic frequency » 0 indicates a characteristic component; when 0 <C |/%fc € | < tt/2, the 

characteristic oscillation frequency is intermediate; arid tt/2 < \0^h^\ < n characterizes a noncharacteristic 
solution component. The last column marked “No” shows the number of target-grid cycles (including those 
two from the FMG algorithm) required to get an approximation with the L x norm of the residual error less 
than 10- 10 . 

Some conclusions derived from analyzing the numerical results are following. 

1* For characteristic components, the target-grid algebraic error is already less than the discretization 
error on Step 2 (obtaining the initial target-grid solution approximation). This is due to the explicit 
terms introduced to all the coarse-grid discretizations. In this wav, we obtain (nearly) the same 
characteristic component discretization errors on all the grids (the difference is proportional to h 4 ). 

2. For components fast oscillating in the cross-characteristic direction, the algebraic error is much better 
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Table 4.1 

Multigrid six-level FMG solver: characteristic arid intermediate solution components. 


t 


e y h 

0 r k 

Discr. 

error 

Algebraic error 

No 

initial 

1 cycle 

2 cycle 

Characteristic components. 

0.1 

0.000 

0.393 

-0.039 

0.09572 

0.01173 

0.00176 

0.00028 

22 

0.3 

0.000 

0.393 

-0.118 

0.11295 

0.03494 

0.00550 

0.00099 

22 

0.5 

0.000 

0.393 

-0.196 

0.04570 

0.06286 

0.00718 

0.00121 

21 

0.7 

0.000 

0.393 

-0.275 

0.11295 

0.08332 

0.01299 

0.00229 

22 

0.9 

0.000 

0.393 

-0.353 

0.09572 

0.10806 

0.01584 

0.00248 

24 

0.1 


1.571 

-0.157 

0.99990 

0.08726 

0.02518 

0.00618 

21 

0.3 


1.571 

-0.471 

1.00006 

0.12234 

0.06842 

0.01987 

24 

0.5 

0.000 

1.571 

-0.785 

0.99996 

0.18733 

0.13551 

0.01981 

23 

0.7 

0.000 

1.571 

-1.100 

1.00006 

0.12716 

0.06039 

0.01920 

23 

0.9 

0.000 

1.571 

-1.414 

0.99990 

0.19453 

0.01811 

0.00537 

21 

Intermediate components. 

0.1 

0.197 

0.393 

0.158 

0.01842 

0.11203 

0.02097 

0.00223 

24 



0.393 

0.087 

0.02952 

0.10400 

0.01697 

0.00325 

22 

0.5 

0.220 


0.023 

0.02903 

0.09289 

0.01502 

0.00428 

20 

0.7 

0.240 


-0.035 

0.02328 

0.14294 

0.02503 

0.00389 

22 

EE1 


0.393 

-0.089 

0.01640 

0.15169 

0.03876 

0.00543 

25 

m 


BUS 

0.040 

0.33494 

0.21031 

0.05228 

0.00902 

21 

EO 


mm 

-0.266 

0.54903 

0.35938 

0.06013 

0.01651 


0.5 

0.220 



0.39790 

0.42716 

0.10210 

0.02764 

22 

0.7 

0.240 

1.571 

-0.860 

0.24156 

0.35623 

0.10455 

0.02441 

23 

0.9 

0.264 

1.571 

-1.150 

0.12979 

0.27564 

0.03267 

0.00633 

22 



0.393 

0.750 

0.14177 

1.07711 

0.26069 

0.07341 

26 





0.16808 

1.43491 

0.33600 

0.08341 

25 





0.20206 

1.67709 

0.66490 

0.09648 

24 




0.684 

0.23211 

2.15520 

0.84395 

0.19261 

26 

0.9 




0.27629 

1.81925 

0.59167 

0.26728 


0.1 

0.789 

mm 

0.632 

0.33086 

0.76440 

0.12903 

0.02681 

24 




0.349 

0.65990 

0.93222 

0.18361 

0.03148 

25 

0.5 

0.878 

1.571 

0.093 

0.60069 

1.09042 

0.24374 

0.04832 

23 

0.7 

0.959 

1.571 

-0.141 

0.51163 

1.04272 

. . 

0.25873 

0.06088 

24 

0.9 

1.057 

1.571 

-0.357 

0.31980 

1.35260 

0.31203 

0.04456 

25 


than the discretization one after the first improving target-grid cycle. 

3. The second improving cycle is needed only for some components moderately oscillating in both the 
characteristic and cross-characteristic directions. On finer grids, these components are eventually 
transferred to the characteristic components. This justifies using FMG algorithm with only one V T 
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Table 4.2 

Multigrid six-level FMG solver: noncharacteristic solution components. 


t 


Oyh 

0 x h 

Discr. 

error 

Algebraic error 

No 

initial 

1 cycle 

2 cycle 

Noncharacteristic components. 

0.1 

1.579 

0.393 

1.539 

0.76817 

2.26395 

0.98437 

0.17369 

28 

0.3 

1.640 

0.393 

1.522 

0.87096 

2.42296 

1.14909 

0.41368 

26 

0.5 

1.756 

0.393 

1.560 

1.04311 

2.91390 

1.88121 

0.64483 

26 

0.7 

1.917 

0.393 

1.643 

1.31669 

3.18658 

1.67748 

0.29945 

27 

0.9 

2.113 

0.393 

1.760 

1.70706 

3.33678 

1.97781 

0.47206 

29 

0.1 

1.579 

1.571 

1.422 

0.85619 

1.85721 

0.38794 

0.07397 

24 

0.3 

1.640 

1.571 

i 

1.169 

1.07279 

2.06405 

0.40828 

0.11672 

25 

0.5 

1.756 

1.571 

0.971 

1.53469 

2.59033 

0.96062 

0.27004 

24 

0.7 

1.917 

1.571 

0.818 

1.59207 

! 2.37166 

0.76083 

0.22509 

25 

0.9 

2.113 

1.571 

0.700 

1.59522 

2.33271 

0.75898 

0.14614 

27 


cycle on coarse levels. 

Note 2: In nonelliptic problems, there are some “pathological 5 * 7 ' noncharacteristic components which 
exhibit very small discretization errors. Noncharacteristic components usually possess relatively large dis- 
cretization errors (compared with characteristic components). However, a very special choice of parameters 
(solution component L and angle of attack <f > ) can result in vanishing discretization errors. It is clear that 
in such special situations we cannot expect the algebraic error to be smaller than or comparable to the 
discretization error at any stage of the algorithm. In spite of the fact that the algorithm fails to reach the 
discretization accuracy for these components, the total (algebraic plus discretization) error in these excep- 
tional cases is much smaller than in neighboring regular cases. Moreover, upon any reasonable perturbation 
the behavior becomes normal: the algebraic error after the two improving target-grid cycles is already sub- 
stantially below the level of the discretization error. It is thus clear in any case that the statement that 
the algebraic error of the FMG solution is less than the discretization error will most likely hold in any real 
calculations (where mostly nonpathological components and angles of attack exist). A detailed analysis of 
this phenomenon (for another type of nonelliptic equations) can be found in [8]. 

5. Adaptive Multigrid Algorithm (AM A). In this section, we present an adaptive multigrid al- 

gorithm approximating the true solution of the differential problem with a relative accuracy e defined in 
advance. The choice of the target grid is a part of the solution process. We restricted ourselves to consider- 
ing only uniform target grids and the target-grid discretization ((2.12) for rn = 1) is also unchanged. Thus, 
the only parameter to be controlled is the target-grid mesh spacing h. The algorithm starts on a very coarse 
grid where the problem is easy to solve and, then, proceeds to finer target grids. On each target grid, the 
FMG algorithm defined in Section 4 is performed to solve the problem. The AMA stops further calculations 
when the relative difference (in a required norm) between the solutions on the two currently finest target 
grids is less than e. In our numerical experiments, the L 0 0 norm of the difference between the solutions on 
the previous target grid and the (injected) current-target-grid solution was served as the stopping criterion. 

Formally, the adaptive multigrid algorithm AMA-FMG employing the FMG cycle described in Section 
4 can be defined in the following 3 steps. 
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FlO . 5.1. AM A — F M G : dependence on solution frequencies 


Step 1 Coarsest- grid solution . The problem on the coarsest grid is solved by some method. 

Step 2 Solution on current target grid. Let u 2h be the solution on the previous target grid with meshsize 
2 h. The current target grid is chosen to be the uniform grid with meshsize h. The FMG algorithm 
is employed to obtain the target-grid solution u h . 

Step 3 Comparison of the solutions. The solution u h is restricted to the previous target grid by some fine- 
to-coarse intergrid transfer operator If/ 1 . In the simplest case, n h is the injection operator. The 
relative difference d r is calculated as 


d'p 


,2h 


j2h„,h\ 


i2h I 


If d r < f, then the solution on the current target grid is considered as final, otherwise algorithm 
proceeds to the next (finer) target grid (Step 2). 

Figures 5.1 and 5.2 demonstrate the algorithm performance for different solution components and dif- 
ferent angles of attack. In these experiments, the known true solution of the differential problem (2. 2), (2.3) 
was U (x, y) = sin (6*1 j- + O^y). The characteristic frequency is calculated as = {9 X + + t' 2 . In 

all the tests shown on Figure 5.1, the nonalignment parameter t was set to t = 0.2 and the frequencies 0$ 
and 0-2 were varied. We tested a large variety of frequencies. Figure 5.1 demonstrates just a representative 
sample of experiments. In the second set of experiments illustrated by Figure 5.2, the frequencies 0$ and 
9-i were fixed but the angle of attack was changed. In all the tests, our goal was to obtain a 1%-accurate 
solution approximation u‘ l (e = 0.01), where the accuracy measured as the L x norm of the relative total error 
(||t/ - a" || 

OO /Ill'll oo < c). Recall, that the algorithm itself does not use the true solution at all. The decision 
whether to stop calculation or proceed to the next finer target grid is made by means of comparison of the 
computed solutions on the current and previous target grids. On the figures, the vertical coordinate marks 
(in the logarithmic scale) the total error of approximations obtained at different stages of the algorithm. 
The vertical lines separate the calculations performed on different, target grids. The coarsest grid problem is 
solved by the FMG algorithm. The first value on each grid is the total error of the approximation obtained 
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^ = 0; 0 2 =16n 



Fig. 5.2. AM A — F M G : dependence on angles of attack 



FlCi. 5.3. AM A — FMG: 2e tolerance 


after the FMG solution interpolation (Step 2 in the FMG algorithm description in Section 4). The two 
following values are the total errors of the approximations obtained after the first and the second target-grid 
cycles respectively. 

The demonstrated AMA performance is quite satisfactory. The total (relative) accuracy of final solutions 
is always much better than the desired 1%-accuracv. This overshooting is actually natural since in order to 
verify that the accuracy of the solution on a grid with meshsize h, the algorithm must compute the solution 


23 




on the next fine grid with meshsize h/2. We tried to change the tolerance of the algorithm weakening the 
stopping criterion to d r < De (D > 1). We tested the same test-cases as in experiments traced on Figure 5.1. 
For D - 2. the algorithm sharply improved its efficiency (see Figure 5.3) finishing the calculations exactly 
on the grids where the required 1%-accuracy had been reached. Further increase of D results sometimes in 
final solutions with too large total errors. 


6. Conclusions and Further Developments. The main difficulty appearing in solving nonelliptic 
equations by multigrid methods is a poor coarse-grid approximation to the characteristic error compo- 
nents. A novel approach proposed in this paper is based on the idea of using semicoarsening together with 
introduction of explicit terms into coarse-grid discretizations in order to maintain the saihe coarse-grid cross- 
characteristic interaction as on the target uniform grid. This constructions allows us to get a good coarse-grid 
approximation to the characteristic components and in this way to solve the aforementioned problem. 

Several multigrid algorithms were tested in this paper. All of them solve the two-dimensional constant- 
coefficient, narrow second-order discretization of the convection equation on uniform Cartesian grids where 
the grid lines do not align with the characteristic direction. The following features of the tested multigrid 
algorithms were reported: 

1 . The algorithms used colored relaxation schemes on all the levels. It makes them very attractive for 
parallel computing. 

2. The residual asymptotic convergence rate of the proposed 1 7 (0,2) multigrid cycle is about 3 per 
cycle. This convergence rate far surpasses the theoretical limit (4/3) predicted for standard multigrid 
algorithms using full coarsening. The reported efficiency does not deteriorate with increasing the 
cycle depth (number of levels) and/or refining the target-grid mesh spacing. 

3. The full multigrid algorithm (FMG) using two V'(0,2) cycles on the target grid and just one V(0,2) 
cycle on all the coarse grids always provides an approximate solution with the algebraic error less 
than the discretization error. The estimates of the total work in the FMG algorithm are between 18 
and 30 minimal work units (depending on the target discretization). Thus, the overall efficiency of 
the FMG solver closely approaches (if does not achieve) the goal of the textbook multigrid efficiency. 

4. A novel approach to deriving a discrete solution approximating the true continuous solution with 
a given relative accuracy is developed. An adaptive multigrid algorithm (AMA) using comparison 
of the solutions on two successive target grids to estimate the accuracy of the current target-grid 
solution is defined. This new criterion for the discrete approximation accuracy is much more effective 
and reliable than the residual monitoring widely used in practice. A desired relative accuracy e 
(0 < f -C 1) is accepted as an input parameter. The final target grid on which this accuracy can be 
achieved is chosen automatically in the solution process. The relative accuracy of the discrete solution 
approximation obtained by AMA is always better than the required e-accuracy. The computational 
work required to compute this approximate solution is (nearly) optimal (comparable with the cost 
of the FMG algorithm applied to solve the problem on the optimally spaced target grid). 

6.1. Extension to Three-Dimensional Problems. In 3D constant-coefficient, case, the characteris- 
tics of the differential equation are still straight lines. Therefore, the low-dimensional prototype is essentially 
the same as (2.4) 


( 6 . 1 ) 


ihy/l + tl+t* 


+3Ui,, ( :, 5u;,_i,j 2 _( y)l 3_ <s + U tl -2,i 2 -2( v ,t 3 -2/..^ =/i 1 ,i 2 ,» 3 , 
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where t y — tan ((f> y ) is the tangent of the angle between the a>axis and the projection of the characteristic 
going through the grid node (ii.hjh) onto the z-y coordinate plane; t z — tan(<^ ; ) is the same for the x-z 
coordinate plane. The “horizontal inclination assumption” is replaced in 3D with the x-axis inclination 
assumption: \t y \ < 1 and \t z \ < 1. Similar to 2D, the full dimensional discretization is derived from (6.1) by 
replacing values at the points with fractional indexes by (bilinear) interpolation from the values at genuine 
grid points placed in the same y-z grid plane. The second-order upwind discretizations are used at the 
outflow boundary. The multigrid solvers employ the semicoarsening procedure, where the only ^-directional 
meshsize is doubled at each coarsening step. Explicit correcting terms introduced on coarse grids include 
discrete approximations to all possible third and fourth derivatives with respect to y and 2. The coefficients 
of these terms are chosen from the condition that the first differential approximation to the coarse-grid 
discretization taken for the characteristic components is exactly the same as the characteristic-component 
FDA to the target-grid discretization. The 3D relaxation schemes used in the multigrid cycles are colored 
plane and pointwise. defect- correction- type schemes where relaxation of a y-z plane (with given ^-coordinate 
i\) replaces the 2D vertical line relaxation wherever it is necessary. The coefficient A of the additional 
stabilizing term which is the discretized fourth derivative with respect to the characteristic variable £ is 
kept to be A = 0.084. Note, that in plane relaxations, the precise solutions of planes are not required. The 
smoothing rate of a 3D plane relaxation scheme employing just one 2D V cycle to solve a plane problem is 
very much the same as the smoothing rate of a relaxation scheme solving plane problems to zero residuals. 
The intergrid transfers used in different multigrid algorithms are characteristic aligned and essentially repeat 
those in 2D: the restriction operators are upwind first order in V r cycles and symmetric second order in FMG; 
the prolongation operators are symmetric second order in V cycles and symmetric fourth order in FMG. 
Preliminary experiments confirm that the efficiency of the proposed approach does not deteriorate in 3D. 

6.2. Extension to Variable Coefficients. A generalization of the presented approach to smooth 
nonrecirculating variable velocity fields can be done in the way first tested in [6]. The cornerstone of this 
technique is a flexible recursive intergrid (fme-to-coarse) parameter transfer providing the target-grid accu- 
racy in tracing the characteristic trajectories on coarse grids. In this way, we can construct an accurate basic 
coarse-grid discretization well aligned with the characteristic track. The cross-characteristic interaction in 
this discretization is again weaker than on the target grid and, therefore, we can supply it with explicit terms 
to get a good coarse-grid correction to the fine-grid characteristic error components. The main changes in 
discretizations to be obtained are coming from the luck of symmetry in the discrete low-dimensional proto- 
type. It is still a one-dimensional four-point second-order discretization of the advection equation but on a 
nonuniform grid. The upwind-biased discretization corresponding to the Van Leer’s scheme with k . = 0 is 
the average of the second order central L° and pure upwind L b schemes: 


(6.2) L c u UM 
L v u UA . 2 


= \{L c u il _ ln + L v u Ut ,\, 


“ hi+hz h x u O + U2 + (* 


i+*i) + (*£ - T^) u h,h - 


— 1 [" / ^ 2+^3 h 2 \ 7 . . + ,, I h. 2 1 . 

~ / l 3 l \ h 2 ^2 + ^3 / ll ’* 2 /*2 1**2 + (*2 + » 2 ) / l 2 +/*3 U *1 _2 * *2 + ( *3 + «3 ) J * 


where h\, /i*2, and h% are distances between the discretization nodes of the low-dimensional prototype mea- 
sured along the characteristic going through the grid point where the discrete operator is defined; k\, 
and £3 are integers denoting the vertical displacements (in meshsizes) from the point (ii,i 2 ); .sq, ,s 2 . and 
S3 (0 < si,s 2 ,s 3 < 1) are the tuning parameters. (See Figure 6.1 for a pictorial explanation of the two- 
dimensional discretization stencil.) Generally speaking, this 2D scheme is second-order accurate only on those 
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® • low-dimcnsionaJ prototype stencil; 0 * 2D discretization stencil. 

Fig. 6.1. Coarse-grid discretization: variable velocity field. 

fine grids where the smoothness of velocity field can be exploit. The coefficients of the explicit terms are now 
calculated from the comparison of characteristic- component FDA coefficients of the true cross-characteristic 
derivatives (rather than the y-directional derivatives as in the constant- coefficient cases) because velocity 
directions may be different for different grid nodes (especially, on coarse grids). Preliminary numerical tests 
demonstrate a high efficiency of the approach in application to the variable-coefficient problems. 

6.3. Avoiding Line Relaxations on Coarse Grids. If the characteristics of the differential equation 
change their general orientation over different parts of the domain, the entire domain should be divided into 
(possibly overlapped) subdomains (each occupying an 0(1) part of the domain and having a unique reference 
axis compatible throughout with the characteristic orientation) and the relaxation sweeps should be applied 
separately on each of the subdomains. In this view, using line relaxation schemes is undesirable because of 
possible luck of a global definition for lines. In many cases, this problem can be easily avoided since there 
is no need to solve simultaneously all the equation centered at the same vertical grid line. In other words, 
a block relaxation updating just a part of the line points at a time would be efficient as well. To keep the 
efficiency, the size of the overlap between neighboring blocks should be proportional to the relative coupling 
RC value (see Sections 3.2.3 and 3.3.1). This is one possible way to adjust the approach to a multiblock 
structure. 

Another, even more efficient way is to widen the basic coarse-grid discretization on the grids (or even at 
separate grid nodes) where the value of the relative coupling is greater than (or comparable with) 1. The 
widening of discretization schemes is illustrated on Figure 6.2. 

W idened discretization schemes possess stronger inherent cross-characteristic interaction and, therefore, 
the weight of the explicit terms in the total cross-characteristic interaction is reduced (RC is smaller). This 
technique is efficient only on grids with large enough RC , since the necessary condition to keep the efficiency 
is that the pointwise relaxation scheme is sensitive to the characteristic error components fast oscillating in 
the cross-characteristic direction. Combination of semicoarsening with widening discretization stencils allows 
us to keep RC bounded (around 1) and avoid in this way using line relaxation schemes making pointwise 
schemes always efficient. 
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Fig. 6.2. Coarse-gnd discretization: widening discretization stencil . 



Fig. 6.3. Residual convergence history in multilevel l 7 {0,2) cycles employing the downstream relaxation. The target grids 
are uniform grids with meshsizes h T = hy = h. The problem is homogeneous (f = 0). uj is the frequency of the incoming 
oscillation . The nonalignment parameter t = 0.2. 


6.4. Note on Downstream Relaxation. The main subject of this paper is multigrid algorithms using 
colored relaxation schemes and, therefore, possessing a great parallelization potential. In this section, we 
remark about V cycles with line downstream (sequential order) relaxation schemes. We already mentioned 
that the efficiency of a colored relaxation scheme improves when more colors (in horizontal direction) are 
used. In this extend, the downstream relaxation is an extreme case where the number of colors coincides 
with the number of grid nodes in the ^-direction. The downstream line relaxation demonstrates the best 
smoothing properties among other schemes. The smoothing factor Snid of this scheme is defined as 
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L{0) - (L b {6 a )e at * + L 4 (6 y )e^) ‘ 

LsWv) + L-i(9,)e- i9 * + L x {9 y )e-™- ' 

See all the definitions (for 9,0 i ,L(§),L j (0 tt ),j = 1,2, 3, 4, 5) in Section 3.3. The absolute value of Sm d (0) is 
always less than 0.4. Figure 6.3 shows excellent residual convergence rate histories of the V(0, 2) cycle using 
the downstream line relaxation scheme. Different plots correspond to cycles starting on different uniform 
target grids. 


S 771(1(0) = max 

jr / 2 <|0el<* 
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error. Estimates of the total work in the FMG algorithm are ranged between 18 and 30 minimal work units 
(depending on the target discretization). Thus, the overall efficiency of the FMG solver closely approaches (if does 
not achieve) the goal of the textbook nmltigrid efficiency. 3) A novel adaptive multigrid approach to deriving a 
discrete solution approximating the true continuous solution with a relative accuracy given in advance is developed. 
The computational complexity of this method is (nearly) optimal (comparable with the complexity of the FMG 
algorithm applied to solve the problem on the optimally spaced target grid). 
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